function [moments paths] = solve_simulate(est_params,glob,param,age_option,wage_option)

    %{
    This function solves and simulates the model for given parameter values
    in est_params. It returns the set of moments we aim to match in the 
    data, as well as the simulation paths.

    How to use input "age_option":
    - set equal to -1 to simulate group whose dist. matches the sample
    - set equal to another number to get a group of people age_sim + 20
    years old

    How to use input "wage_option":
    - set equal to 0 to keep wage fixed across simulation
    - set equal to 1 for PERMANENT change in wage (for Mashallian
    elasticity)
    - set equal to 2 for TRANSITORY change in wage (for Frisch elasticity)

    %}

    param.beta     = est_params(1);        % Consumption weight
    param.delta    = est_params(2);        % Discount factor
    %param.gammah   = est_params(3);        % Hours limit
    %est_params

    %% Initialize: set up for backward induction

    % matrix to store coefficients/values/policies in
    c       = zeros(glob.Ns,glob.T-1);
    vals    = zeros(glob.Ns,glob.T-1);
    pols    = zeros(glob.Ns,glob.T-1);
    cons    = zeros(glob.Ns,glob.T-1);
    hrs     = zeros(glob.Ns,glob.T-1);
    bc      = zeros(glob.Ns,glob.T-1);

    % compute utility at terminal state
    uT = param.beta*log((1-glob.tau)*glob.at-glob.gammac) + (1-param.beta)*log(glob.gammah);

    % compute first set of basis coefficients by hand (asset choice will be 0)
    RHS = param.beta*log(glob.s+(1-glob.tau)*glob.at-glob.gammac) + (1-param.beta)*log(glob.gammah);
    c(:,glob.T-1) = funbas(glob.fspace,glob.s)\RHS;

    %% Solve backwards

    for t=glob.T-2:-1:1
        obj = @(ap)valfunc(ap,glob.s,c(:,t+1),t,param,glob);

        % choose bounds for optimization
        if t>= glob.R
            ub = (1+glob.r).*(glob.s+(1-glob.tau)*glob.at);
        else
            ub = (1+glob.r).*(glob.s+(1-glob.tau)*glob.wt(t).*glob.gammah);
        end

        Ap = goldenx(obj,glob.amin.*ones(glob.Ns,1),ub);
        pols(:,t) = Ap;
        cons(:,t) = menufun('consumption',Ap,glob.s,t,param,glob);
        hrs(:,t) = menufun('hours',Ap,glob.s,t,param,glob);
        vopt = obj(Ap);
        vals(:,t) = vopt;
        c(:,t) = funbas(glob.fspace,glob.s)\vopt;         % compute coefficients

        % check that the budget constraint binds
        if t>= glob.R
            bc(:,t) = Ap - (1+glob.r).*(glob.s + (1-glob.tau)*glob.at - cons(:,t));
        else
            bc(:,t) = Ap - (1+glob.r).*(glob.s + (1-glob.tau)*glob.wt(t).*hrs(:,t) - cons(:,t));
        end

    end


    %% Simulate a panel: lottery winners

    % set up
    rng(222);
    wealth_nowin = nan(glob.N,glob.T);
    wealth_win   = nan(glob.N,glob.T);
    hours_nowin = nan(glob.N,glob.T);
    hours_win   = nan(glob.N,glob.T);
    earn_nowin = nan(glob.N,glob.T);
    earn_win   = nan(glob.N,glob.T);
    age          = nan(glob.N,glob.T);

    if age_option == -1
        % draw initial ages
        % approximate histogram from the paper:
        lb = 1;
        for a=1:length(glob.age_bins)
            ub = lb + glob.age_probs(a)*glob.Na;
            age(lb:ub,1) = ones(ub-lb+1,1).*glob.age_bins(a);
            lb = ub;
        end
        age = age(1:end-1,:);
    else
        age(:,1) = age_option;
    end

    % set up pre-lottery wealth distribution: everyone at the mean
    wealth_nowin(:,1) = ones(size(wealth_nowin(:,1))).*glob.mean_wealth;

    % draw lottery prizes
    wealth_win(:,1)   = wealth_nowin(:,1) + randsample(glob.prizes,glob.N,true,glob.prize_dist)';
    dA0               = wealth_win(:,1) - wealth_nowin(:,1);

    % now we're ready to simulate
    %tic
    for t = 2:glob.T

        % update the age
        age(:,t) = age(:,t-1) + 1;

        % update wealth conditional on winning and not winning, for those alive
        for a = 1:glob.T-1

            % find the indices of people at that age
            inds = find(age(:,t)==a);
            if isempty(inds)
                continue
            end
            lowbound = min(inds);
            uppbound = max(inds);

            % take their wealth levels: points at which to evaluate policies
            w_nowin = wealth_nowin(lowbound:uppbound,t-1);
            w_win   = wealth_win(lowbound:uppbound,t-1);
            wpoints = [w_nowin; w_win];
            
            % policy interpolant
            pol_a         = funfitxy(glob.fspace,glob.s,pols(:,a)); 
            w_t1          = funeval([pol_a],glob.fspace,[wpoints]);

            % compute policy: solve exactly
%             obj = @(ap)valfunc(ap,wpoints,c(:,a),a,param,glob);
%             if a >= glob.R
%                 ub = (1+glob.r).*(wpoints+glob.at);
%             else
%                 ub = (1+glob.r).*(wpoints+glob.wt(a).*param.gammah);
%             end
%             w_t1 = goldenx(obj,glob.amin,ub);

            % update tomorrow's wealth
            wealth_nowin(lowbound:uppbound,t) = w_t1(1:end/2);
            wealth_win(lowbound:uppbound,t)   = w_t1(end/2+1:end);

            % compute hours and earnings
            hours_nowin(lowbound:uppbound,t) = menufun('hours',wealth_nowin(lowbound:uppbound,t),...
                wealth_nowin(lowbound:uppbound,t-1),a,param,glob);
            hours_win(lowbound:uppbound,t) = menufun('hours',wealth_win(lowbound:uppbound,t),...
                wealth_win(lowbound:uppbound,t-1),a,param,glob);
            earn_nowin(lowbound:uppbound,t) = menufun('earnings',wealth_nowin(lowbound:uppbound,t),...
                wealth_nowin(lowbound:uppbound,t-1),a,param,glob);
            earn_win(lowbound:uppbound,t) = menufun('earnings',wealth_win(lowbound:uppbound,t),...
                wealth_win(lowbound:uppbound,t-1),a,param,glob);

        end

    end
    %toc
    
    %% Simulate a panel: wage changes
%     if wage_option ~= 0
%         
%         % set up
%         rng(222);
%         wealth_wc       = nan(glob.N,glob.T);
%         hours_wc        = nan(glob.N,glob.T);
%         earn_wc         = nan(glob.N,glob.T);
%         age             = nan(glob.N,glob.T);
% 
%         if age_option == -1
%             % draw initial ages
%             % approximate histogram from the paper:
%             lb = 1;
%             for a=1:length(glob.age_bins)
%                 ub = lb + glob.age_probs(a)*glob.Na;
%                 age(lb:ub,1) = ones(ub-lb+1,1).*glob.age_bins(a);
%                 lb = ub;
%             end
%             age = age(1:end-1,:);
%         else
%             age(:,1) = age_option;
%         end
%         
%         % set up wealth distribution
%         wealth_wc(:,1)   = ones(size(wealth_wc(:,1))).*glob.mean_wealth;
%         
%         % change the wages
%         inds = (0:length(glob.wt)-1)';
%         if wage_option == 1 & age_option == -1
%             indices = (inds >= 1);
%         elseif wage_option == 2 & age_option == -1
%             indices = (inds == 1);
%         elseif wage_option == 1 & age_option ~= -1
%             indices = (inds >= age_option + 1);
%         elseif wage_option == 2 & age_option ~= -1
%             indices = (inds == age_option + 1);
%         end
%         glob.wt = glob.wt + indices.*glob.wt*0.1;
%         
%         % now we're ready to simulate
%         %tic
%         for t = 2:glob.T
% 
%             % update the age
%             age(:,t) = age(:,t-1) + 1;
% 
%             % update wealth conditional on winning and not winning, for those alive
%             for a = 1:glob.T-1
% 
%                 % find the indices of people at that age
%                 inds = find(age(:,t)==a);
%                 if isempty(inds)
%                     continue
%                 end
%                 lowbound = min(inds);
%                 uppbound = max(inds);
% 
%                 % take their wealth levels: points at which to evaluate policies
%                 w_wc   = wealth_wc(lowbound:uppbound,t-1);
%                 wpoints = [w_wc];
% 
%                 % policy interpolant
%                 pol_a         = funfitxy(glob.fspace,glob.s,pols(:,a)); 
%                 w_t1          = funeval([pol_a],glob.fspace,[wpoints]);
% 
%                 % compute policy: solve exactly
%     %             obj = @(ap)valfunc(ap,wpoints,c(:,a),a,param,glob);
%     %             if a >= glob.R
%     %                 ub = (1+glob.r).*(wpoints+glob.at);
%     %             else
%     %                 ub = (1+glob.r).*(wpoints+glob.wt(a).*param.gammah);
%     %             end
%     %             w_t1 = goldenx(obj,glob.amin,ub);
% 
%                 % update tomorrow's wealth
%                 wealth_wc(lowbound:uppbound,t)   = w_t1;
% 
%                 % compute hours and earnings
%                 hours_wc(lowbound:uppbound,t) = menufun('hours',wealth_wc(lowbound:uppbound,t),...
%                     wealth_wc(lowbound:uppbound,t-1),a,param,glob);
%                 earn_wc(lowbound:uppbound,t) = menufun('earnings',wealth_wc(lowbound:uppbound,t),...
%                     wealth_wc(lowbound:uppbound,t-1),a,param,glob);
% 
%             end
% 
%         end
%         
%         paths.wealth_wc = wealth_wc;
%         paths.hours_wc  = hours_wc;
%         paths.earn_wc   = earn_wc;
%         paths.wages     = glob.wt;
%         
%     end
% 

    %% compute moments

    % effects for pre-tax earnings
    diff_earnings = earn_win - earn_nowin;
    avg_eff_earnings = zeros(10,1);
    for i = 1:length(avg_eff_earnings)
        eff_earnings = (diff_earnings(:,i+1))./dA0;
        avg_eff_earnings(i) = nanmean(eff_earnings)*100;
    end

    % average hours conditional on winning
    %Hbar = nanmean(hours_win(:,2:11),2);

    % final moments
    %moments = [avg_eff_earnings; nanmean(Hbar)];
    moments = avg_eff_earnings;

    % pack up simulation paths
    paths.age           = age;
    paths.wealth_nowin  = wealth_nowin;
    paths.wealth_win    = wealth_win;
    paths.hours_nowin   = hours_nowin;
    paths.hours_win     = hours_win;
    paths.earn_nowin    = earn_nowin;
    paths.earn_win      = earn_win;

end